# recuperation de l'ensemble des donnees
create_summary_data <- function(inbfile)
{
	# création d’un fichier de déséquilibre de Hardy Weinberg (hwe.log, hwe.hh, hwe.hwe)
	system(paste("plink --noweb --bfile ",inbfile," --hardy --out hwe",sep=""))

	# Création d'un fichier de différence de call rate (% de génotypage) entre cas et contrôles (dmiss.log, dmiss.hh, dmiss.missing)
	system(paste("plink --noweb --bfile ",inbfile," --test-missing --out dmiss",sep=""))

	#Création d'un fichier de fréquence allélique (maf.log, maf.hh, maf.frq)
	system(paste("plink --noweb --bfile ",inbfile," --freq --out maf",sep=""))
	
	# Mise en commun des résultats 
	hwe <- read.table(file="hwe.hwe",h=T)
	dmiss <- read.table(file="dmiss.missing",h=T)
	maf <-  read.table(file="maf.frq",h=T)

	data <- merge(hwe,dmiss, by=c("CHR","SNP"))
	names(data)[9] <- "P.HW"
	names(data)[12] <- "P.dmiss"
	data <-merge(data,maf, by=c("CHR", "SNP", "A1", "A2"))
	write.table(data, file="summary_data.txt",quote=F,sep="\t")
}

#creation d'une liste de SNPs à éliminer
snp_selection <- function(HW, miss, maf)
{
	# chargement des donnees
	data <- read.table(file="summary_data.txt", h=T)
	aff <- data[data$TEST=="AFF",]
	unaff <- data[data$TEST=="UNAFF",]
	tmp <- merge(aff,unaff,by=c("CHR","SNP"))
	# Elimination des SNPs HW < HW
	tmp_aff <- tmp[tmp$P.HW.x < HW & !is.na(tmp$P.HW.x),"SNP"] 
	tmp_unaff <- tmp[tmp$P.HW.y < HW & !is.na(tmp$P.HW.y),"SNP"]

	# Elimination des SNPs avec un % de missing  > miss
	tmp_aff2 <- tmp[tmp$F_MISS_A.x > miss & !is.na(tmp$F_MISS_A.x),"SNP"]
	tmp_unaff2 <- tmp[tmp$F_MISS_U.y > miss & !is.na(tmp$F_MISS_U.y),"SNP"]

	# Elimination des SNPs avec MAF < maf
	tmp_aff3 <- tmp[tmp$MAF.x < maf & !is.na(tmp$MAF.x),"SNP"]
	tmp_unaff3 <- tmp[tmp$MAF.y < maf & !is.na(tmp$MAF.y),"SNP"]

	# fusion des donnees
	allSNP <- c(as.character(tmp_aff),as.character(tmp_aff2),as.character(tmp_aff3),as.character(tmp_unaff),as.character(tmp_unaff2),as.character(tmp_unaff3))
	allSNP <- unique(allSNP)

	# Enregistrement des données
	write.table(allSNP,file="exclusionSNP.txt",quote=F, row.names=F)
}

# elimination des SNPs – creation de nouveau fichiers plink
exclusionSNP <- function(inbfile)
{
	# pour les cas
	system(paste("plink --noweb --bfile ",inbfile,"  --exclude exclusionSNP.txt --filter-cases --make-bed --out cases",sep=""))
	# pour les controles
	system(paste("plink --noweb --bfile  ",inbfile," --exclude exclusionSNP.txt --filter-controls --make-bed --out controls",sep=""))
}

# creation des fichiers cluster - a faire pour les cas et les controles
create_cluster_file <- function(plate_file,  status) #status : 1 for cases, 2 for controls
{
	# chargement et préparation des données
	plates <- read.table(file=plate_file,h=T)
	if(status == 1){
	info <- read.table(file="cases.fam")
	}
	if(status == 2){
	info <- read.table(file="controls.fam")
	}
	info <- info[-c(3,4,5,6)]
	names(info)[1] <- "FID"
	names(info)[2] <- "ADN"

	# fusion des donnees
	cluster <- merge(info,plates, by=c("ADN"))
	cluster <- cluster[-c(3)]
	names(cluster)[1] <- "IID"
	cluster<- cluster[c(2,1,3)]

	# enregistrement des données
	if(status == 1) {
	write.table(cluster,file="cluster-cases.dat",quote=F, row.names=F)
	}
	if(status == 2){
	write.table(cluster,file="cluster-ctrls.dat",quote=F, row.names=F)
	}
}


# creation du fichier de phenotype - a faire pour les cas et les controles
create_pheno_file <- function(cluster_file, status) #status : 1 for cases, 2 for controls
{
	# chargement des données
	cluster <- read.table(file=cluster_file,h=T)
	plate_list <-sort(unique(cluster$NPLAK))
	x <- length(cluster$FID)
	y <- length(plate_list)

	# création de la matrice des données
	mat <- matrix(1,nrow=x,ncol=y)
	nplak <- cluster$NPLAK

	# remplissage de la matrice de données
	for (i in 1:x){
	j=1
	repeat{
	if(nplak[i]==plate_list[j])break else{j=j+1;}
	}
	mat[i,j] <- 2
	}

	# mise en forme de la matrice
	pheno <- cbind(cluster[-c(3)],mat)

	# nommer les colonnes
	names <- colnames(cluster[-c(3)])
	col_names <- c(names,plate_list)
	colnames(pheno) <- col_names
	num <- colnames(pheno[-c(1,2)])

	# enregistrement des données
	if(status == 1) {
	write.table(pheno,file="pheno_cases.txt",quote=F, row.names=F, col.names=T)
	write.table(num,file="header-pheno-file1.txt")
	}
	if(status == 2){
	write.table(pheno,file="pheno_ctrls.txt",quote=F, row.names=F, col.names=T)
	write.table(num,file="header-pheno-file2.txt")
	}
}


# creation des fichiers d'association sous plink
association <- function()
{
	# pour les cas
 	system(paste("plink --noweb --bfile cases --pheno pheno_cases.txt --all-pheno --assoc --out cases_plate",sep=""))

	# pour les contrôles
	system(paste("plink --noweb --bfile controls --pheno pheno_ctrls.txt --all-pheno --assoc --out ctrls_plate",sep=""))
}


info_chisq_pval <- function(p.names="P", status=NULL){ #status : 1 for cases, 2 for controls
	
	if(status == 1){
		#chargement du fichier de phenotype
		assoc_name <- "cases_plate"
		num <- read.table(file="header-pheno-file1.txt")
	}
	if(status == 2){
		#chargement du fichier de phenotype
		assoc_name <- "ctrls_plate"
		num <- read.table(file="header-pheno-file2.txt")
	}	

	# chargement des donnees de plaque
	for ( i in 1:length(p.names)) {
		name <- paste(assoc_name, num[i,1], sep=".")
		name <- read.table(file=paste(name,"assoc", sep="."),h=T)
		if(i == 1) { 
			pval <- name[c(1,2)]
			chisq <- name[c(1,2)]
		}
		tmp_pval <- name$P
		tmp_chisq <- name$CHISQ

		#  fusion des données
		pval <- cbind(pval,tmp_pval)
		chisq <- cbind(chisq,tmp_chisq)
	}

	#  on nomme les colonnes
	col_names <- c("CHR","SNP",p.names)
	colnames(pval) <- col_names
	colnames(chisq) <- col_names
	
	# moyenne des chisq par plaque 
	chisq_tmp <- chisq[-c(1,2)]
	mean_chisq<- as.vector(apply(X=chisq_tmp,na.rm=T , MARGIN=2, FUN=mean))

	# minimum des p-value par plaque
	pval_tmp <- pval[-c(1,2)]
	min_pval<- as.vector(apply(X=pval_tmp,na.rm=T , MARGIN=2, FUN=min))

	# creation d'un fichier info_chisq_pval
	info_chisq_pval <- cbind(p.names,mean_chisq,min_pval)	
	write.table(info_chisq_pval, file="info_chisq_pval.txt", quote=F, row.names=F)

	# ajout de la p-value minimum
	P_min <- apply(X=pval_tmp, MARGIN=1, FUN=min)
	pval <- cbind(pval,P_min)

	# enregistrement des données
	if(status == 1){
		write.table(pval,"p-values-cases.txt", quote=F, row.names=F)
	}
	if(status == 2){
		write.table(pval,"p-values-controls.txt", quote=F, row.names=F)
	}
}


# creation d'un QQplot par plaque - a faire pour les cas et pour les controles
qqplot.p<- function (p.names="P", status=NULL) #status : 1 for cases, 2 for controls
{
	# chargement de la fonction plot.QQs
	source("qqmultiple_function.R")

	#  visualisation des p-values des plaques 
	if(status == 1){
		pval <- read.table(file="p-values-cases.txt",h=T)
		png("qqplot-plaque-cases.png", width=940, height=940, res=72)
	}
	if(status == 2){
		pval <- read.table(file="p-values-controls.txt",h=T)
		png("qqplot-plaque-controls.png", width=940, height=940, res=72)
	}

	nb_plate <- length(p.names)
	nb_plot <- ifelse(nb_plate%%4 ==0, length(p.names)%/%4, (length(p.names)%/%4)+1)
	x<-1
	y<-4
	par(mfrow=c(2,2))
	
	for (j in 1:nb_plot){
		plot.QQs(data1=pval,pobs.name=p.names[x:y])
		x<-x+4
		if(y+4 <= length(p.names)){y <- y+4}
		else {y<- length(p.names)}
	}
	dev.off()
}

# creation d'un QQplot par plaque - a faire pour les cas et pour les controles
qqplot.p<- function (p.names="P", status=NULL) #status : 1 for cases, 2 for controls
{
	# chargement de la fonction plot.QQs
	source("qqmultiple_function.R")
	
	if(status == 1){
		#chargement du fichier de phenotype
		assoc_name <- "cases_plate"
		num <- read.table(file="header-pheno-file1.txt")
	}
	if(status == 2){
		#chargement du fichier de phenotype
		assoc_name <- "ctrls_plate"
		num <- read.table(file="header-pheno-file2.txt")
	}	

	# chargement des donnees de plaque
	for ( i in 1:length(p.names)) {
		name <- paste(assoc_name, num[i,1], sep=".")
		name <- read.table(file=paste(name,"assoc", sep="."),h=T)
		if(i == 1) { pval <- name[c(1,2)]}
		tmp <- name$P
		#  fusion des données
		pval <- cbind(pval,name$P)
	}

	#  on nomme les colonnes
	col_names <- c("CHR","SNP",p.names)
	colnames(pval) <- col_names
	
	# ajout de la p-value minimum
	pval_tmp <- pval[-c(1,2)]
	P_min <- apply(X=pval_tmp, MARGIN=1, FUN=min)
	pval <- cbind(pval,P_min)

	# enregistrement des données
	if(status == 1){
		write.table(pval,"p-values-cases.txt", quote=F, row.names=F)
	}
	if(status == 2){
		write.table(pval,"p-values-controls.txt", quote=F, row.names=F)
	}

	#  visualisation des p-values des plaques 
	if(status == 1){
		png("qqplot-plaque-cases.png", width=940, height=940, res=72)
	}
	else {
		png("qqplot-plaque-controls.png", width=940, height=940, res=72)
	}

	nb_plate <- length(p.names)
	nb_plot <- ifelse(nb_plate%%4 ==0, length(p.names)%/%4, (length(p.names)%/%4)+1)
	x<-1
	y<-4
	par(mfrow=c(2,2))
	
	for (j in 1:nb_plot){
		plot.QQs(data1=pval,pobs.name=p.names[x:y])
		x<-x+4
		if(y+4 <= length(p.names)){y <- y+4}
		else {y<- length(p.names)}
	}
	dev.off()
}

## creation d'un manhattan plot pour chaque plaque
#mhtplot.p <- function (status=NULL) #status : 1 for cases, 2 for controls
#{
#	if(status == 1){
#		pval <- read.table(file="p-values-cases.txt", h=T)
#	}
#	if(status == 2){
#		pval <- read.table(file="p-values-controls.txt", h=T)
#	}
#	colors <- c("red","blue","green","cyan","yellow","gray","magenta","red","blue","green",          "cyan","yellow","gray","magenta","red","blue","green","cyan","yellow","gray","magenta","red")
#	 mhtplot(data,control=mht.control(colors=colors))
#}

# mesure de l'excès de valeurs significatives
lambdaGC <- function (p.names="P", status){ #status : 1 for cases, 2 for controls
	
	if(status == 1){
		#chargement des noms des colonnes
		assoc_name <- "cases_plate"
		num <- read.table(file="header-pheno-file1.txt")
	}
	if(status == 2){
		#chargement des noms des colonnes
		assoc_name <- "ctrls_plate"
		num <- read.table(file="header-pheno-file2.txt")
	}	
	# chargement des donnees de plaque
	for ( i in 1:length(p.names)) {
		name <- paste(assoc_name, num[i,1], sep=".")
		name <- read.table(file=paste(name,"assoc", sep="."),h=T)
		if(i == 1) {chisq <- name$CHISQ}
		if(i>1){
		tmp_chisq <- name$CHISQ
		#  fusion des données
		chisq <- cbind(chisq,tmp_chisq)
		}
	}

	#  on nomme les colonnes
	colnames(chisq) <- p.names

	# moyenne des chisq par plaque 
	mean_chisq<- apply(X=chisq,na.rm=T , MARGIN=2, FUN=mean)

	# enregistrement des donnees
	if(status == 1){
	write.table(as.data.frame(mean_chisq), file="mean_chisq_cases.txt", quote=F, row.names=F)
	}
	if(status == 2){
	write.table(as.data.frame(mean_chisq), file="mean_chisq_controls.txt", quote=F, row.names=F)
	}
}

#test <- function(tab){
#for(row in 1:length(tab[1])){
#for (i in 3:13) {
#if(tab[row,i] == tab[row,14]){
#name_min[row] <- names(tab)[i]
#i = i+1 
#}
#}
#}
#}

#test <- function(tab){
#	Plak_min <- rep(0,length(tab$CHR))
#	for (i in 1:length(tab$CHR)){
#		j=3
#		repeat{
#			if(tab[i,j]== tab[i,14] & j <14 )break else{j <- j+1}
#		}
#	Plak_min <- names[j]
#	}
#}

#plak <- rep("NA",length(cases$CHR))
#for(i in 1:length(cases$CHR)){ 
#j=3
#print(i)
#repeat{
#if(j<14 & cases[i,j]==cases[i,14]) break else{j <- j+1}
#}
#plak[i] <- names(cases)[j]
#}


# mesure de l'excès de valeurs significatives
lambdaGC <- function (p.names="P", status){ #status : 1 for cases, 2 for controls
	
	if(status == 1){
		#chargement des noms des colonnes
		assoc_name <- "cases_plate"
		num <- read.table(file="header-pheno-file1.txt")
	}
	if(status == 2){
		#chargement des noms des colonnes
		assoc_name <- "ctrls_plate"
		num <- read.table(file="header-pheno-file2.txt")
	}	
	# chargement des donnees de plaque
	for ( i in 1:length(p.names)) {
		name <- paste(assoc_name, num[i,1], sep=".")
		name <- read.table(file=paste(name,"assoc", sep="."),h=T)
		if(i == 1) {chisq <- name$CHISQ}
		if(i>1){
		tmp <- name$CHISQ
		#  fusion des données
		chisq <- cbind(chisq,tmp)
		}
	}

	#  on nomme les colonnes
	colnames(chisq) <- p.names

	# moyenne des chisq par plaque 
	mean_chisq<- apply(X=chisq,na.rm=T , MARGIN=2, FUN=mean)

	# enregistrement des donnees
	if(status == 1){
	write.table(as.data.frame(mean_chisq), file="mean_chisq_cases.txt", quote=F, row.names=F)
	}
	if(status == 2){
	write.table(as.data.frame(mean_chisq), file="mean_chisq_controls.txt", quote=F, row.names=F)
	}
}

#source("analyse_plaque.R")
#create_summary_data("racont_final")
#snp_selection(0.0001,0.05,0.1)
#exclusionSNP("racont_final")
#create_cluster_file("../Plaks_N.txt",1)
#create_cluster_file("../Plaks_N.txt",2)
#create_pheno_file("cluster-cases.dat",1)
#create_pheno_file("cluster-ctrls.dat",2)
#association()
#info_chisq_pval(p.names=c("P9cases","P10","P11","P12","P13","P14","P15","P16","P17","P18","P19"), status=1)
#info_chisq_pval(p.names=c("P1","P2","P3","P4","P5","P6","P7","P8","P9ct"), status=2)
#qqplot.p(p.names=c("P9cases","P10","P11","P12","P13","P14","P15","P16","P17","P18","P19"),status=1)
#qqplot.p(p.names=c("P1","P2","P3","P4","P5","P6","P7","P8","P9ct"),status=2)

#lambdaGC(p.names=c("P9cases","P10","P11","P12","P13","P14","P15","P16","P17","P18","P19"),status=1)
#lambdaGC(p.names=c("P1","P2","P3","P4","P5","P6","P7","P8","P9ct"),status=2)
